#include <stdio.h>
#include <math.h>

//编译：g++ -o run_fit_curve_v2 fit_curve_v2.c -lm
//编译运行：g++ -o bin/run_fit_curve_v2 c/fit_curve_v2.c -lm && ./bin/run_fit_curve_v2
//对比第一版，只是模块化更好了而已，性能没有变

// 高斯消元法函数
void gaussianElimination(double matrix[3][4], int n) {
    double factor;
    int i, j, k;
    // 高斯消元
    for (k = 0; k < n - 1; k++) {
        // 寻找主元
        int max_row = k;
        for (i = k + 1; i < n; i++) {
            if (fabs(matrix[i][k]) > fabs(matrix[max_row][k])) {
                max_row = i;
            }
        }
        if (max_row != k) {
            double temp[4];
            for (j = 0; j <= n; j++) {
                temp[j] = matrix[k][j];
                matrix[k][j] = matrix[max_row][j];
                matrix[max_row][j] = temp[j];
            }
        }

        for (i = k + 1; i < n; i++) {
            factor = matrix[i][k] / matrix[k][k];
            for (j = k; j <= n; j++) {
                matrix[i][j] -= factor * matrix[k][j];
            }
        }
    }

    // 回代求解
    for (i = n - 1; i >= 0; i--) {
        for (j = i + 1; j < n; j++) {
            matrix[i][n] -= matrix[j][n] * matrix[i][j];
        }
        matrix[i][n] /= matrix[i][i];
    }
}
// 函数来求解抛物线方程的系数
void solveParabola(double points[3][2], double *coefficients) {
    double matrix[3][4];

    // 构建增广矩阵
    for (int i = 0; i < 3; i++) {
        matrix[i][0] = points[i][0] * points[i][0];
        matrix[i][1] = points[i][0];
        matrix[i][2] = 1;
        matrix[i][3] = points[i][1];
    }

    // 调用高斯消元法函数
    gaussianElimination(matrix, 3);

    // 提取结果
    coefficients[0] = matrix[0][3] * 100000;
    coefficients[1] = matrix[1][3] * 100000;
    coefficients[2] = matrix[2][3] * 100000;
	printf(">>> 计算参数：%f %f %f \r\n",coefficients[0],coefficients[1],coefficients[2]);
}


int main(){

    //拟合曲线使用
    double points[3][2] = {
        {35658,9.808}, // 点1
        {67369,15.67}, // 点2
        {78092,19.261999}  // 点3
    };
    double coefficients[3]={0,0,0};//方程a、b、c值,y=x

    //处理函数
    solveParabola(points, coefficients);

    printf(">>> y = %lfx^2 + %lfx + %lf\n", coefficients[0], coefficients[1], coefficients[2]);
    float x = 166387.0;
    printf("x=%f so y=%f \r\n",x,((coefficients[0]*pow(x,2))+(coefficients[1]*x)+coefficients[2])/100000);


    return 0;
}


